﻿using System;

namespace Plugins.Samples
{
    /// <summary> 
    /// Calculation of PI using the BBP (Bailey-Borwein-Plouffe) formula. 
    /// </summary> 
    public static class Pi
    {
        /// <summary> 
        /// Powers of two. Used in the binary exponentiation algorithm (Knuth Vol II - 4.6.3). 
        /// </summary> 
        static int[] s_powsOf2;

        /// <summary> 
        /// Hexadecimal characters. 
        /// </summary> 
        static string HEX = "0123456789ABCDEF";

        /// <summary> 
        /// Static constructor. 
        /// </summary> 
        static Pi()
        {
            s_powsOf2 = new int[25];

            // 
            // Set powers of two for use in binary exponentiation algorithm. 
            // 
            s_powsOf2[0] = 1;
            for (int i = 1; i < s_powsOf2.Length; i++)
                s_powsOf2[i] = s_powsOf2[i - 1] * 2;
        }

        /// <summary> 
        /// This method calculates the fractional part of SUM[i = 0..+INF](1 / (16^i * (8 * i + j))). 
        /// </summary> 
        /// <param name="j">See formula.</param> 
        /// <param name="pos">Starting position for requested hexadecimal part of PI's fraction.</param> 
        /// <returns>See formula.</returns> 
        public static double Sum(int j, int pos)
        {
            /* 
             * Implementation based on the following observation: 
             * 
             *   SUM[i = 0..+INF](1 / (16^i * (8 * i + j))) 
             * 
             *   =   SUM[k = 0..d]((16^(d - k) mod (8 * k + j)) / (8 * k + j))) 
             *     + SUM[k = d + 1..+INF](16^(d - k) / (8 * k + j))) 
             */

            double res = 0.0;

            int k = 0;

            // 
            // First part of the sum using binary exponentiation with modulo operation. 
            // 
            for (; k < pos; k++)
            {
                double n = 8 * k + j;
                res += PowMod(16.0, pos - k, n) / n;
                res -= (int)res;
            }

            // 
            // For the remainder terms, go till we reach the precision of 10^-17, 
            // which allows to reach the required precision in hexadecimals. 
            // 
            for (; ; k++)
            {
                double t = Math.Pow(16.0, pos - k) / (8 * k + j);
                if (t < 1e-17)
                    break;
                res += t;
                res -= (int)res;
            }

            return res;
        }

        /// <summary> 
        /// Returns a specified number raised to the specified power, module the given modulus. 
        /// </summary> 
        /// <param name="x">A double-precision floating-point number to be raised to a power.</param> 
        /// <param name="y">A double-precision floating-point number that specifies a power.</param> 
        /// <param name="m">A double-precision floating-point number that specified the modulus.</param> 
        /// <returns>x^y mod m</returns> 
        static double PowMod(double x, double y, double m)
        {
            if (m == 1.0)
                return 0.0;

            // 
            // We look for the largest bit position used by the exponent. 
            // E.g. 16^9 mod 9 = (16 * 16^(2^3)) mod 9 
            //                               ^ 
            // 
            int i;
            for (i = 0; i < s_powsOf2.Length; i++)
                if (s_powsOf2[i] > y)
                    break;

            double pow2 = s_powsOf2[i - 1];

            // 
            // Set remaining number of exponentiations and result to initial values. 
            // 
            double rem = y;
            double res = 1.0;

            for (int j = 0; j < i; j++)
            {
                // 
                // "Remainder" steps that do not fit in binary exponentiation. 
                // 
                // j   res             res' 
                // --  --------------  ---- 
                // 0   16^0 mod 9 = 1  16^1 mod 9 = (1 * 16) mod 9 = 7 
                // 3   16^8 mod 9 = 4  16^9 mod P = (4 * 16) mod 9 = 1 
                // 
                if (rem >= pow2)
                {
                    res = (res * x) % m;
                    rem -= pow2;
                }

                pow2 /= 2;

                // 
                // Binary exponentation step. 
                // 
                // j   res             res' 
                // --  --------------  ------------------------------ 
                // 0   16^1 mod 9 = 7  16^2 mod 9 = (7 * 7) mod 9 = 4 
                // 1   16^2 mod 9 = 4  16^4 mod 9 = (4 * 4) mod 9 = 7 
                // 2   16^4 mod 9 = 7  16^8 mod 9 = (7 * 7) mod 9 = 4 
                // 
                if (pow2 >= 1.0)
                {
                    res = (res * res) % m;
                }
            }

            return res;
        }

        /// <summary> 
        /// Converts the fractional part of the specified double-precision floating-point number 
        /// to hexadecimal representation. 
        /// </summary> 
        /// <param name="d">Double-precision floating-point number in base 10.</param> 
        /// <returns>String representation of the number's fractional part in base 16.</returns> 
        /// <example> 
        /// Sample: d := 0.14159265358979334 
        /// 
        ///    f := (d - (int)d )   16 * f               Hex((int)(16 * f)) 
        ///    -------------------  -------------------  ------------------ 
        /// 0. 0.14159265358979334   2.2654824574366934  2 
        /// 1. 0.26548245743669340   4.2477193189870945  4 
        /// 2. 0.24771931898709450   3.9635091037935126  3 
        /// 3. 0.96350910379351260  15.4161456606962020  F 
        /// 4. 0.41614566069620200   6.6583305711392313  6 
        /// 5. 0.65833057113923130  10.5332891382277010  A 
        /// 6. 0.53328913822770100   8.5326262116432190  8 
        /// 7. 0.53262621164321900   8.5220193862915039  8 
        /// 8. 0.52201938629150390   8.3523101806640625  8 
        /// 9. 0.35231018066406250   5.6369628906250000  5 
        /// 
        /// Reverse: h := 0.243F6A8885 
        /// 
        ///    Weight       Value 
        ///    -----------  -----------  
        /// 0. 2 * 16^-1    0.12500000000000000000 
        /// 1. 4 * 16^-2    0.01562500000000000000 
        /// 2. 3 * 16^-3    0.00073242187500000000 
        /// 3. F * 16^-4    0.00022888183593750000 
        /// 4. 6 * 16^-5    0.00000572204589843750 
        /// 5. A * 16^-6    0.00000059604644775391 
        /// 6. 8 * 16^-7    0.00000002980232238770 
        /// 7. 8 * 16^-8    0.00000000186264514923 
        /// 8. 8 * 16^-9    0.00000000011641532183 
        /// 9. 5 * 16^-10   0.00000000000454747351 
        ///                +---------------------- 
        ///                 0.14159265358921400000 
        ///                   ^^^^^^^^^^^^ 
        /// </example> 
        public static string ToHex(double d)
        {
            char[] c = new char[10];

            for (int i = 0; i < c.Length; i++)
            {
                d = 16 * (d - (int)d);
                c[i] = HEX[(int)d];
            }

            return new string(c);
        }
    }
}
